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1 Introduction 



Qh ' Photoacoustic Imaging is one of the recent hybrid imaging techniques, which 

attempts to visualize the distribution of the electromagnetic absorption coef- 
ficient inside a biological object. In photoacoustic experiments, the medium 
is exposed to a short pulse of a relatively low frequency electromagnetic (EM) 
wave. The exposed medium absorbs a fraction of the EM energy, heats up, 
and reacts with thermoelastic expansion. This induces acoustic waves, which 
can be recorded outside the object and used to determine the electromag- 
netic absorption coefficient. The combination of EM and ultrasound waves 
^ ' (which explains the usage of the term hybrid) allows one to combine high 

contrast in the EM absorption coefficient with high resolution of ultrasound. 
The method has demonstrated great potential for biomedical applications, 
including functional brain imaging of animals [WPKXS03] , soft-tissue charac- 
terization, and early stage cancer diagnostics [KKMROO], as well as imaging 
of vasculature [ZLB07]. For a general survey on biomedical applications see 
' [XW06]. In comparison with the X-Ray CT, photoacoustics is non- ionizing. Its 

further advantage is that soft biological tissues display high contrasts in their 
ability to absorb frequency electromagnetic waves. For instance, for radiation 
in the near infrared domain, as produced by a Nd:YAG laser, the absorp- 
' tion coefficient in human soft tissues varies in the range of 0.1/cm— 0.5/cm 

[CPW90]. The contrast is also known to be high between healthy and can- 
cerous cells, which makes photoacoustics a promising early cancer detection 
technique. Another application arises in biology: Multispectral optoacoustic 
tomography technique is capable of high- resolution visualization of fluorescent 
proteins deep within highly light-scattering living organisms [RDVMPKN]. 
In contrast, the current fluorescence microscopy techniques are limited to the 
depth of several hundred micrometers, due to intense light scattering. 

Different terms are often used to indicate different excitation sources: Op- 
toacoustics refers to illumination in the visible light spectrum, Photoacoustics 
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is associated with excitations in the visible and infrared range, and Ther- 
moacoustics corresponds to excitations in the microwave or radio-frequency 
range. In fact, the carrier frequency of the illuminating pulse is varying, which 
is usually not taken into account in mathematical modeling. Since the corre- 
sponding mathematical models are equivalent, in the mathematics literature, 
the terms opto-, photo-, and thermoacoustics are used interchangeably. In this 
article, we are addressing only the photoacoustic tomographic technique PAT 
(which is mathematically equivalent to the thermoacoustic tomography TAT). 

Various kinds of photoacoustic imaging techniques have been implemented. 
One should distinguish between photoacoustic microscopy (PAM) and tomog- 
raphy (PAT). In microscopy, the object is scanned pixel by pixel (or voxel by 
voxel). The measured pressure data provides an image of the electromagnetic 
absorption coefficient [ZMSW06]. Tomography, on the other hand, measures 
pressure waves with detectors surrounding completely or partially the object. 
Then the internal distribution of the absorption coefficients is reconstructed 
using mathematical inversion techniques (see the sections below). 

The common underlying mathematical equation of PAT is the wave equa- 
tion for the pressure 



dj_ { Mabs(x)/3(x)J(x) 
dt { ' V c p (x) 
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Here c p denotes the specific heat capacity, J is the spatial intensity distribu- 
tion, ^i a b s denotes the absorption coefficient, denotes the thermal expansion 
coefficient and Co denotes the speed of sound, which is commonly assumed 
to be constant. The assumption that there is no acoustic pressure before the 
object is illuminated at time t — is expressed by 



p(x,t) = 0, x€l 3 ,i<0 



(2) 



In PAT, j(t) approximates a pulse, and can be considered as a <5-impulse 8{t). 
Introducing the shorthand notations 



p(x) := 



//abs(x)/3(x)J(x) 

c p (x) 



(3) 



one reduces (1) and (2) to 



1 d 2 p 



(x,i) - V 2 p(x,i) = 0, xeM 3 ,t>0. 



(4) 



with initial values 



p(x,0)=p(x), ^(x,0)=0 xel 3 



(5) 
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The quantity p in (1) and (3) is a combination of several physical parameters. 
All along this paper p should not be confused with the source term 



In PAT, some data about the pressure p(x, t) are measured and the main task 
is to reconstruct the initial pressure p from these data. While the excitation 
principle is always as described above and thus (4) holds, the specific type of 
data measured depends on the type of transducers used, and thus influences 
the mathematical model. 

Nowadays there is a trend to incorporate more and more modeling into 
photoacoustic. In particular, taking into account locally varying wave speed 
and attenuation. Even more there is a novel trend to qualitative photoacoustics, 
which is concerned with estimating physical parameters from the imaging 
parameter of standard photoacoustics. In this paper we focus on attenuation 
correction, where we survey some recent progress. Inversion with varying wave 
speed has been considered for instance in [AK07, HKN08], and is not further 
discussed here. 

The outline of this paper is as follows: First, we review existing attenuation 
models and discuss their causality properties, which we believe to be essential 
for algorithms for inversion with attenuated data. Then, we survey causality 
properties of common attenuation models. We also derive integro-differential 
equations which the attenuated waves are satisfying. In addition we discuss 
the ill-conditionness of the inverse problem for calculating the unattenuated 
wave from the attenuated one. 



2 Attenuation 

The difficult issue of effects of and corrections for the attenuation of acoustic 
waves in PAT has been studied [RZA06, BGHNP07, PG06, KSBIO], although 
no complete conclusion on the feasibility of these models has been reached. 

Mathematical models for describing attenuation are formulated in the fre- 
quency domain, taking into account that attenuation disperses high frequency 
components more rapidly over traveled distance. Let Q (x, t) denote the atten- 
uated wave which originates from an impulse ((^-distribution) at x = at 
time t = 0. In mathematical terms Q is the Green- function of attenuated wave 
equation. Moreover, we denote by 



/(x,i) = -£(t)p(x) 



x e M 3 ,t > o . 



(6) 




the Green function of the unattenuated wave equation; that is, it is the solution 
of (4), (5) with constant sound speed c(x) = c and initial conditions 
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e?o(x,0) = and ^( x ,0)=£ x>t . 
Common mathematical formulations of attenuation assume that 



(x,w) = cxp (-/3*(|x| ,w)) J"{£o} (x,w) , x G R 3 , w G 



(8) 



Here T {•} denotes the Fourier transform with respect to time t (cf. Ap- 
pendix 8). Applying the inverse Fourier transform J 7-1 {•} to (8) gives 



G(x-,t) = K(x.,t) *t <?o( x , t) (*t time convolution) 



where 



?~ X {exp (-/3*(|x| , •))} (t) . 



(9) 
(10) 



From (9) and (7) it follows that 

S(x,i) =A"(x,t) * t £o(x,t) 

= / K(x,t-r) - , ? rfr 



4-7T Ixl 



A" 



47T Ixl 



Consequently, 



e?(x, i + |x| /c ) = A(x, t)/(47T |x|) . 



(11) 



Moreover, we emphasize that the Fourier transform of a real and even (real 
and odd) function is real and even (imaginary and odd). Since Q and Go are 
real valued, K must be real valued and consequently the real part 9?(/3*) of (3* 
has to be even with respect to the frequency to and 3(/3*) has to be odd with 
respect to w. Attenuation is caused if 3ft(/3*) is positive and since then j3* has 
a nonzero imaginary part due to the Kramers-Kronig relation, attenuation 
causes dispersion. In the literature the following product ansatz is commonly 
used 

(12) 



<«>) = «» M 



u) G R, x G 



In the sequel we concentrate on these models and use the following terminol- 
ogy: 

Definition 1. We call (3* of standard form if (12) holds. Then the function 

a*:R^C (13) 



is called standard attenuation coefficient and a = s R(a* ) is called the attenu- 
ation law. We also call (3* the attenuation coefficient. 
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From the relation (12), it follows that SR(a*) is even, 3(a*) is odd, and 3?(a*) > 
(the last inequality guarantees attenuation). 

In the following we summarize common attenuation coefficients and laws: 
In what follows ao denotes a positive parameter and 



«o = 



cos (§7) 

is a possibly non-positive coefficient. 



(0 < 7 g N) 



(14) 



• Frequency Power Laws: 

Let < 7 ^ N. The frequency power law attenuation coefficient is 
defined by 

a * P i( u ) = «o(-iw) 7 = a M 7 (cos (| 7 ) - isgn(w) sin (7^7)) (15) 
for uj e R. Therefore, the attenuation law is given by 

(16) 



a p i(uj) = a \lo\ 



These models have been considered for instance in [S94, S95, WHBM00, 
WMM05] . 

Let 7 = 1 and ujq ^ 0, the attenuation coefficient is defined by 



a*;(w) := a \u\ +i-a w log 



7T 



Ul 

UJQ 



The attenuation law is 



a p i(ui) := a M . 



(17) 



(18) 



This model has been considered in [S95, WMM05]. 
• Szabo: Let < 7 ^ N. The attenuation coefficient 1 of Szabo's law is 
defined by 



w) = — v^-iw) 2 + 2a co(-iw)T +1 + i 
co 



co 



(19) 



We denote Szabo's attenuation law by 



a sz (oj) := SR(a* z (o;)) 



For small frequencies a sz (w) behaves like ao |w| 7 . This model has been 
considered in [S94, S95] where, in addition, also a model for 7 e N has 
been introduced. 



1 In this paper the root of a complex number is always the one with non-negative 
real part. 
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Thermo- Viscous Attenuation Law: (see e.g. [KFCSOO, S94]): Here, 
for r > 0, the attenuation coefficient is defined by 



c4H = 



— 1W 1U 

c VI — i^o w c o 



with attenuation law 




(20) 



(21) 



For small frequencies c<t„(w) behaves like : rr i£ -. That is the thermo- viscous 



law approximates a power attenuation law with exponent 2. 
• Nachman, Smith and Waag [NSW90]: Consider a homogeneous and 
isotropic fluid with density po in which N relaxation processes take place. 
Then the attenuation coefficient of the model in [NSW90] reads as follows: 



CO 



1 N 1 

Co 1_ l~lT m M _ 

M Jv ^ 1 - i7VnW 



(22) 



All parameters appearing in (22) are positive and real. K m and r m denote 
the compression modulus and the relaxation time of the m— th relaxation 
process, respectively, and 



c := 



c 



and f m ■= r rn (1 - N c% po K m ) • (23) 



' 1 +T lm =l C PO K ™ 

The last two definitions imply that 

S2 



N _ 



Co _ 1 \ ~> T m 

~r2 ~ AT ' 



(24) 



m— 1 



We denote the according attenuation law by 



a nsw (u) := $t(a„ sw (u)) 



• Greenleaf and Patch [PG06] consider for 7 e {1,2} the attenuation 
coefficient 



which, since it is real, equals the attenuation law 



a gp {uj) = ft(a* (w)) 



(25) 



In [NSW90] they use the notion for c and c for cq. 
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• Chen and Holm [CH04]: This model describes the attenuation as a 
function of the absolute value of the vector- valued wave number k E M 3 
(instead of the frequency wgR), Let T ZD denote the 3D— Fourier trans- 
form 

•F 3D {/(k)}(x) = -=2= f exp(ix.k) /(k)dk, 

v( 27r ) J 

then the Green function of the attenuated equation is defined by 

„, . H(t)cl ( / a / \ \ sin(.B(-) t) 1 , . 

g(x, t) = F 3D jcxp (A(-) t) [ B ^> > | (x) (26) 

where, for given ai > 0, 

A(k) := -ai co |k| 7 , B(k) := c ^/|k| 2 - a\ |k| 2 ^ . (27) 

• In [KSB10] we proposed 

a*k.M = F= f ( 7 e(l,2],r >0), (28) 

where the square root is again the complex root with positive real part. 
Let 7 G (1)2]. Then, for small frequencies we have 



a sin(f (7- 1)) 

afcshM w r w| 7 > . 

2 c T 

Thus our model behaves like a power law for small frequencies. 

Distinctive features of unattenuated wave propagation (,i.e. the solution of 
the standard wave equation) are causality and finite wave front velocity. It is 
reasonable to assume that the attenuated wave satisfies the same distinctive 
properties as well. In the following we analyze causality properties of the 
standard attenuation models. 



3 Causality 

In the following we present some abstract definitions and basic notations. In 
the remainder x will always denote a vector in three dimensional space. When 
we speak about functions, we always mean generalized functions, such as for 
instance distributions or tempered distributions - we recall the definitions of 
(tempered) distribution in the course of the paper. 

Definition 2. A function f := /(x, t) defined on the Euclidean space over 
time (i.e. in R 4 ) is said to be causal if it satisfies /(x,t) = for t < 0. 
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Notation & Terminology 1 Let A : D — > D be a linear operator, where 
^ D is an appropriate set of functions from R 4 to ffi. In this paper we 
always assume that A satisfies the following properties: 

• A is shift invariant in space and time. That is, for every function f and 
every shift L := L(x,t) := (x — xo,t — to), with xo € K 3 and to e K, it 
ftoZds t/iat 

^(/ o L) = (Af) o L . 

• A is rotation invariant in space. That is, for every function f and every 
rotation matrix R, it holds that 

A(Rf) = R(Af) . 

• A is causal. That is, it maps causal functions to causal functions. From 
(29) it follows that A is causal, if and only if the associated Green function 
is causal. 

Definition 3. The Green function of A is defined by 

G := g(x,t) = A5 x , t (x,t) . 

Remark 1. The operator A is uniquely determined by Q and vice versa. This 
follows from the fact that 




Af(x , t ) =A \ I J /(x - x, t - to)5 x ,t(x, i) dxdt 

R R3 

/(xq — x, t — to)<?(x, t) dxdt . 



(29) 



R R3 

Moreover, we use the following terminology and abbreviations: 
From the rotation invariance of A it follows that 

f(x) := sup{t : Q(x, r) = for all t < t} , (30) 

is rotationally symmetric, which allows us to use the shorthand notation 

T(r) = f(x) where r = |x| . (31) 

With this notation (30) can be equivalently expressed as 

G {x,t + T(|x|)) = for every t < . (32) 

In physical terms T(|x|) denotes the travel time of a wave front originating 
at position at time t — and traveling to x. 



Photoacoustics 9 

Because Q is rotationally symmetric we can write 

g(x,T(|x|)) = G(r,T(r)) with r = |x| . 

Taking the inverse function of T , which we denote by S = S(t), we then 
find 

g(x,T(\ x \)) = g(S(t),t), 

The wave front is the set 

W := {(x,T(|x|)) :xeR 3 } 

The wave front speed V is the variation of the location of the wave front 
as a function of time. That is, 



T'(r) 



(33) 

r=S(t) 



.Here T' denotes the derivative with respect to the radial component r. 

• VFe say i/iai .A has a finite speed of propagation if there exists a constant 
c such that 

< (T'(r)r 1 < c < oo . (34) 
In this case it follows from (33) that the wave front velocity satisfies 

V(t) < c < oo . (35) 

• We call an operator A strongly causal, if it is causal and satisfies the finite 
propagation speed property. 

The following lemma addresses the case of attenuation coefficients of stan- 
dard form and gives examples of strongly causal operators A. 

Lemma 1. Let /3*(|x| ,cj) = a*(u) |x| be of the standard form (12) and A 
(29) be the operator defined by the Green function Q, which is defined in (9). 
Then A is strongly causal if and only if for every xel 3 the function 

t^_LF-i{ exp (-a»|x|)} , 

V27T 



defined in (10), is causal. 

Proof. We assume that A is strongly causal. It follows from [KSB10, Theorem 
3.1] that there exists a constant c, which is smaller than or equal to the 
wave speed co from (4), which satisfies T(|x|) = ^ for all x e M 3 . Using 
the definitions of the travel time T(|a;|) and (10), it follows from (11) that 
t — > K(x,t) is causal. 

Now, for every x e M 3 let K be causal. Then from (11) it follows that 
t ^ Q (x,t + |x| /co)) is causal. Since T(|x|) denotes the largest positive time 
period for which t — >• Q (x, t + T(|x|)) is causal, we have for r > 0: 
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< - < T(r) = / -J- . (36) 
co 7 V(s) 

o 

Here V is parameterized with respect to the distance s at time t of the wave 
front from its origin. As shown in the proof of [KSB10, Theorem 3.1], the 
fact that f3* is of standard form together with (36) implies that there exist a 
constant c such that T(r) = r/c for all r > 0. But then from (36) it follows 
< r/c < r/c < oo and consequently < c < c < oo. 

Finally we explain the above notation for the standard wave equation: 

Remark 2. In the case of the standard wave equation the wave front is the 
support of the Green function Q , the wave front velocity is c , and T(|x|) = ^ 
denotes the travel time of the wave front. 



4 Strong Causality of Attenuation Laws 

In this section we analyze causality properties of attenuation laws. We split the 
section into two parts, where the first concerns numerical studies to determine 
the kernel function K, defined in (10), and the second part contains analytical 
investigations. 

In Figures 1, 2 and 3 we represent the attenuation kernels according to 
power, Szabo's, and the thermo-viscous law. The figures already indicate 
that power laws with index greater than 1 violate causality. In the following 
we support these computational studies by analytical considerations. Thereby 
we make use of distribution theory, which we recall first. Generally speaking 
Distributions are generalized functions: 

Definition 4. We use the abbreviations: 

• V := Cq°(R, C) is the space of infinitely often differentiable functions from 
R to C which have compact support. 

• S is the space of infinitely often differentiable functions from R to C which 
are rapidly decreasing. A function f : R — > C is rapidly decreasing if for 
all i,j e N 



for \x\ — > oo . 



S is a locally convex space (see [Y95] for a definition) with the topology 
induced by the family of semi-norms 

Pp J (/)=su P P(x)f^(x) 

where P is a polynomial and j € No. The topology on a locally convex set 
is defined as follows: U C S is open, if for every f € U there exists e > 
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Fig. 1. Simulation of K(jx.,t) for the frequency power law with (7,00) G 
{(0.5,0.1581), (1.5,0.0316), (2.7,0.0071), (3.3,0.0027)}, c = 1 and |x| = \. In the 
first example 7 < 1 and thus the function is causal. For all other cases it is non 
causal. 

and a finite non-empty set J' of polynomials and a finite set of indices K' 
such that 

p| {geS:p Pjk (g-f)<e}CU. 

• The space of tempered distributions, S' , is the space of linear continuous 
functionals and S . 

• A functional L : S — > C is continuous if there exists a constant C > and 
a seminorm ppj such that 

\Lu\ < Cpp } j(u) , for every u e S 

(see [Y95, Sect. 1.6, Thm. 1] 

In the following we give some examples of tempered distributions and 
review some of their properties. The examples are taken from [Y95, Sec6.2, 
Ex. 3] and [DL02c, Remark 6]. 

Example 1. (Examples of Tempered Distributions) 

• Let 1 < p < 00 and / G L P (M.,C), then the linear operator T<fi = 
J R f(x)(j)(x)dx is a tempered distribution. In the following we identify / 
and T, and this clarifies the terminology / G S' later on. 





Fig. 2. Simulation of K(x,t) for Szabo's frequency law with (7,00) 6 
{(0.5,0.1581), (1.5,0.0316), (2.7,0.0071), (3.3,0.0027)}, c = 1 and |x| = \. 




Fig. 3. Left: if(x, t) defined by the complex thermo-viscous attenuation law with 
tq = 10~ 5 , Co = 1 and fixed |x| = |. Right: The proposed law (28) for 7 = 2 with 
Qi = 1, To = 10~ 5 , co = 1 and fixed |x| = | is causal. 



• S C. S' - thereby already the above relation between functions and tem- 
pered distributions is used. 

• Distributions with compact support are tempered distributions. For in- 
stance the ^-Distribution is a tempered distribution. 

• Polynomials are tempered distributions. 
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• The functions / of L ; 1 oc (R) which are uniformly bounded by a polynomial 
for \x\ sufficiently large, are tempered distributions. 3 

Lemma 2. • The pointwise limit f : R — > C of a sequence of functions 
{f n : R — >• C} C S' , is again a tempered distribution. 

• Let f e S', thenF{f} G 5' and J" -1 {/} G 5'. 

In the following we review Theorem 4 on p294 ff from [DL02b] which 
characterized when a generalized function / G <S'(R) is causal, that is, when 
SU PP(/) [0, oo). Below we use the following notation 

C e := {z G C : > e} . 

Theorem 2. (Theorem 4 on p294 ff in [DL02b]) Let f G S'(R). Then f is 
causal 4 if and only if 

1. There exists a function F : Co := {£ + i 77 : 77 > 0} — > C, which is 
holomorphic in the interior Co := {£ + i 77 : 77 > 0}. 5 

U. for a/Z /iired 77 > and (el, F(£ + i 77) is a tempered distribution with 
respect to the variable £ and /or 77 — > F(( + i 77) is convergent (with 
respect to the weak topology on S' ) to T {/} : R — > C. 6 

3. For every e > 0, t/iere exists a polynomial P such that 

\F(z)\ < P(\z\) for ZGC £ . 

Remark 3. The definition of the Fourier transform in this chapter has a differ- 
ent sign as in [DL02b] and consequently also Co denotes the upper half plane 
and not the lower half plane as in [DL02b] . 

For analyzing attenuation laws, we use the following corollary, which is 
derived from Theorem 2. 

Corollary 1. Let a* : R — > C be continuous and let there exist a holomorphic 
extension to Co, which for the sake of simplicity of notation is again denoted 
be a* . In addition, let a*(£ + irj) — > for 77 — > pointwise. We denote 

by a : Co — > C the real part of a* . 7 Moreover, we assume that there exists a 
constant C such that 

a(u) >C for allcu G R . (37) 

3 A function / is an element of L ; 1 oc (E) if it is in L 1 on every compact set. 

4 In Theorem 4 on p294 ff in [DL02b] the assumption that / is strongly causal is 
expressed by / G £>+, which is the set of distributions with support in [0, +00). 

5 A function F : Co — > C is holomorphic in Co if it is complex differentiable in Co. 
Sometimes the functions are also refered to as analytic or regular functions or 
conformal maps. 

6 A function F : Co — > C which satisfies Items 1,2 of Theorem 2 is called holomor- 
phic extension of J-{f}. 

7 a : Co — > C extends the function cj G R — > ol(lo) but is not an holomorphic 
extension. 
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1. If in addition 

a(z) > C for all zeC . (38) 
Then, for every x€l 3 , the function 

t -+ K(x, t) := -^-F" 1 {cxp (-«*(•) |x|)} (t) (39) 
V 27T 

is causal. 

2. On the other hand, if there exists C\ > 0, \i > and C 2 £ K and a 
sequence {z n } in Co such that 

a(z n )<-C 1 \z n \»-C 2 , (40) 

then K violates causality. 

Proof. Let x e R 3 fixed. We apply Theorem 2 to /(•) = X(x, •)• Therefore, 
we have 

V^F{/}(w)=exp(-a*(w)|x|) . 

Under the assumption (37), taking into account that a* ; is continuous, 
is in Lj oc (M.) and bounded by a constant polynomial, thus in S 1 (cf. Example 
2), and consequently, according to Lemma 2, / G 5'. Therefore, the general 
assumption of Theorem 2 is satisfied. 

The function z € Co i-> -F(z) := cxp (— a* (2) |x|) is an extension of 
27rJ r {/} (w), which is holomorphic in Co- Thus Item 1 of Theorem 2 holds. 

• For proving the first assertion, it follows from (38) that for z = £ + i?7 with 
77 >0, 

|F(*)|<exp(-f7|x|), (41) 

which, in particular, shows that for all rj > 0, the functions £ — > F(£ + i 77) 
is a tempered distribution (cf. Lemma 2). Hence Item 1 of Theorem 2 
holds. 

Moreover, since by assumption a*(£ + in) — > for 77 — > pointwise, 

F(£ + ir/) converges to F(£) — 27rJ r {/} (£) pointwise. Because the limit is 
a tempered distribution and the convergence is with respect to the weak 
topology S' (cf. Lemma 2). Hence Item 2 of Theorem 2 holds. 
Moreover, from (41) it follows that |-F(£)| is bounded by a constant poly- 
nomial. Hence Item (3) of Theorem 2 holds and therefore Theorem 2 guar- 
antees that 1 1-> K (x, t) is causal. 

• For the second case, Item 3 of Theorem 2 is violated. Consequently, K(x, •) 
is not causal. 

Power Laws 

Theorem 3. Let < 7 € R, &e not an odd number, and uj G R h-> 
a* z (w) = <S (— icj) 7 fee £/ie power law attenuation coefficient from (15), with 
<So = a o/ cos (§7) as ^ f-^- T/ien, the function K, defined in (39), is causal 
if and only if '7 € (0, 1). 
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Proof. Let x e R 3 be fixed. The function z e C ^ = a (-i2) 7 is the 

holomorphic extension of to e M i-> We prove or disprove causality by 

using Corollary 1. 

For z = \z\ exp (i</>) it follows from (15) that 

K((-iz) 7 ) =%t{\A 1 exp (i 7 (^-tt/2))) = |z| 7 cos( 7 (<?i> - tt/2)) . 

This implies that 

a pl (z) = a 3?((-iz) 7 ) = S cos( 7 (0 - tt/2)) |z| 7 . (42) 

In particular, if z = w € R, then (/> is either or ir. Taking into account the 
definition of «o and that the cos-function is symmetric around the origin, it 
follows that 

Opi(w) = a H 7 >0. (43) 

Thus (37) holds. 

• Let 7 e (0,1): Every z = \z\ exp(i^) e Co satisfies <p G [0, 7r]. Conse- 
quently 7(0 — 7r/2) G [— 7r/2,7r/2] and thus cos(7(0 — 7r/2)) is uniformly 
non-negative. Even more for 7 e (0, 1) the coefficient So, defined in (14), 
is positive. In summary, we have that there exists a constant C\ > such 
that 

a p i(z) > d |z| 7 > for z e C . (44) 

Thus (38) holds and application of Corollary 1 shows that K is causal. 

• Let 7 e (1, 3) U (5, 7) U • • • . Then a < 0. The sequence 

{z„ := nexp(i7r/2) = in}„ eN (45) 
consists of elements of Co and satisfies assumption (40), that is, 

a p i{z n ) = Jx^\z n \' y . (46) 
<o 

Application of Corollary 1 shows that K is not causal. 

• Let 7 e (3, 5) U (7, 9) U • • • . We fix some < 6 < tt/2, and define 

\ 7/2 7 

The sequence 

{z n := n exp(i0)} (47) 

consists of elements of C . Under the above assumptions, it follows that 
(So > and therefore 

Oipl{z n ) = j5q^ cos(7r/ 2 + £) |z„| 7 . (48) 
>o ~~ <o 

Thus form Corollary 1 the assertion follows. 
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In the following we analyze the following family of variants of power laws: 
"pi+M = "of-iw) 1 + ai(-iw) , (49) 
which have been considered in [S95, WHBMOO]. 

Theorem 4. Let < 7 ^ N and a* pl+ as defined in (49). Moreover, let K be 
as in (39). Then, if 7 > 1, K is not causal. For 7 G (0, 1) K is causal if and 
only if ai G [0, 00). 

Proof. The holomorphic extension of lo G M — » a* ;+ (u;) is the function 

a pZ+<» = ao(-iz)'' + ai (-iz) 

and consequently 

a pl+ (z) = a |*f (l + g Izl 1 ^) cos ( 7 (^-|)). 

• For 7 > 1 we have that 1 + §7 |z| 1-7 — > 1 for |z| — > 00. Let {z n } as defined 
in (45) or (47). Then, since for both sequences \z n \ — > 00, it follows from 
(46), (48) that the according sequences {z„} satisfy (40) for n sufficiently 
large, respectively. Thus K is not causal. 

• For 7 G (0, 1) and a\ > the assertion follows already from the fact 
that a p i+(uo) > a p i{uj) and that the later already satisfies (38). Thus K is 
causal. 

• Let 7 G (0, 1) and a\ < 0. Then for some < S < —a\ fixed, we can find 
a constant C2 such that for all z G C 

&o M 7 + ai |z| < (ai +5) |z| - C 2 . 
<o 

Consequently, for {z n = in}, we have 

a p i+(z n ) < (ai + <5) |n| - C 2 . 
which shows (40). Thus K is not causal. 

Power law with 7 = 1 

Theorem 5. Let a pl be as in defined in (17). Then the function K, defined 
in (39) is not causal. 

Proof. First, we prove that 

_ , . . 2ao . { z\ 
zet|)H a i(z) := a Q z + i z log — 

7T \Wo/ 
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is the holomorphic extension of u — > This assertion follows from the 

facts 

<s&M = q4M for w>o, 



lim i — log [ — — — ) = i — log 

J7->0+ 7T \ U>Q J 7T 



w 



- 2 for uj < 



Since 

a p ;(u;) = 5R(a*,(w)) = a M > 0, 

Corollary 1 is applicable. For the elements of the sequence {z n := in}„ G N in 
C 

%l( z n) = "log — 

is real and therefore equals a p i(z n ) and thus (40) holds. Thus Corollary 1 
gives the assertion. 



Szabo's Model: 



Proposition 1. For ao > let a* sz be the coefficient of Szabo's model (19). 
Then, for 7 G (0, 1), the function K (10) is a causal function and for 7 > 1 
with 7 ^ N, K violates causality. 

Proof. Without loss of generality we assume that c = 1. The holomorphic 
extension of a* sz : M -> C from (19) is 



z e C -> «: z (z) - (-iz) [v/l + 25 (-iz)T-i-l 

First, we make some general manipulations which can be used in several 
ways: Let z = £ + i 77 G Co- We use the polar representation 

z=\z\ exp (i 0) , := <f>(z) G [0, 7r] . 

Then 

a: z («) = Nexp(i(0-7r/2))!?(z), 

where 

S 

With this notation we have 

= l + 2a cos((5)|z| 7 " 1 
^(z))=2a sin((5)|zr 1 , 



= l + 2a kr _1 exp(i ( 5) , 
= <5(z) := (0- 7r /2)( 7 -l) 



(50) 



(51) 



(1 + 2<5 cos(<5)) 2 + A& 2 sin((5) 2 
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Representing # in polar coordinates, 

&{z) = &(z) cxp(i 6>(z)) 

we get 



$(z) 



<P{z) exp(i0(z)/2) 



with 



(52) 



6{z) = arctan(9f(!?(«))/R(!?(«)))) G {-it, tt) . 



Note, that \J^{z) is the complex root with non-negative real part, which 
meets the general assumption of the paper. Moreover, we have 

a sz {z) = M(a* sz (z)) = r)M(*(z)) + ■ (53) 

First, we prove that $l(\P(z)) > 0: We use the elementary inequality 

cos(0(z)) < cos 2 (6»(z)/2) , 

and cos(0(z)/2) > which imply that 



R V*(*) 



#(z) cos(0(z)/2) 



#(z) cos 2 (0(z)/2) 



(54) 



> V»(*(*)) 



Now, let 7 e (0, 1). Since $l(&(z)) > 1 for 7 G (0, 1), it follows that for all 

z G Co 



J7&(#(z)) = 775ft ^ ^(z)J - 77 > ri\J^{z)) - 77 > . 
Thus 77&(#(2)) > 0. 

Now we show that z — > ^(^(z)) is uniformly bounded from below by 
in Co- Thus according to (53) a sz is uniformly bounded from below, and 
thus from Corollary 1, it follows that t — > K(x,t) is causal. 
Using the definition of 0, (52), and the facts that 6 G [0, (1 — 7)71-/2] for 
e [0,tt/2] and 5 G [(7 - 1)tt/2,0) for G (7r/2,7r] it follows from the 
monotonicity of tan on (— 7r,7r) that 



0(; 



arctan 



2a s'm(S) \z 



7-1 



1 1 + 2a cos(5) |z| 7 y 
Now, noting that sgn(£) = sgn (sin(0(z)/2)) it follows that 



[0,(5] for all <p G [0,tt/2] , 
[5,0) for all 4> G (tt/2,7t] 



C3(f(z)) = £y ^(2) sin(0(z)/2) > 0. 
Thus the assertion follows from Corollary 1. 



(55) 
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Assume that 7 > 1. Let z = £+i 77 with r\ — 0. Since the square root in (50) 
is such that $l(\P(z)) > for z € Co, property (102) in the Appendix 
implies £3?(<f(z)) > for z = £ and hence a sz (z = £) > 0. Thus (37) holds 
and we can apply Corollary 1. 

Let 76 (1, 3)U(5, 7)U - • • , which implies that cos (777/2) < 0, and con- 
sequently, ao < 0. For sufficiently large n the elements of the sequence 
{z n :— in} satisfy 



a sz (z n ) = ndi (yi - 2 \a Q \ n^ 1 - l) 



< — n . 



which shows that (40) holds with fi = 1, Ci = 1/2 and C2 = 0, and 
hence the assertion follows from Corollary 1. 

Let 7 e (3, 5) U (7, 9) U • • • , which implies that cos (771-/2) > 0, and 
consequently, q.q > 0. Now, let z n := ncxp(i(/>) with 



7T 
+ - . 



7-12 

Since 7 > 3, we have <j) & ( 7r /2,7r), and therefor 3?(z„) = n cos(</>) < 
and Ss(z n ) = n sin(0) > 0. Moreover, 

£(*„) = 1 -2 |a | n^" 1 

and thus for sufficiently large n we have 3?(tf'(z„)) = — 1 and 



3? &(z n )j > Ci n^- 1 ^ 2 for some constant C\ > 0. Hence it follows 
that 

a sz (z n ) = -$S{z n ) + $t(z n ) %(\f<P(z n )) 

= -n |sin(0)| - n |cos(0)| C x n^' 1 ^ 2 
and therefore (40) holds. Thus from Corollary 1 the assertion follows. 

Thermo- Viscous Attenuation Law 

Theorem 6. Let c ,t > and let a* tv as defined in (20). Then the kernel 
function K violates causality. 

Proof. Since the function z e Co — > 1 — iro z does not vanish, the function 

z e Co -> a t *„(z) - ,7 1Z . + — 

coV1-ito^ c 

is the holomorphic extension of u> E M — > 0^(0;). That (37) holds follows 
from the identity (21). For the sequence {z n } ne ^ := {injneN we get for n 
sufficiently large 



* 1 \ n 

Co 



1 

- 1 



VI + t n 



1 

< -7T— n 
~ 2 c 



Thus (40) holds. 

The second part of Corollary 1 implies that K is not causal. 
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Model of Nachman, Smith and Waag 
Theorem 7. Let a^ sw as in (22). If 

f m < r m for all m e {1, . . . , N} , (56) 
then the kernel function K is causal. 

Proof. Since for all z € C and all 1 < m < N, 1 — ir m z does not vanish, 



z e C -s- < aw (2) = 



— 1 z 
co 



\ 



_j_ \ - 1 - 1 T m Z 
N ^ I -IT, 



1 



is the holomorphic extension of u e M — > 0*^(0;). 
We use a similar notation as in Proposition 1. 



z = £ + in = \z\ exp(i^) e C , 



with some <f> e [0,7r]. 



<C™(*) = ^ exp (i (0 - tt/2)) - 1) , 



where 



N 



= \J2 with ■= 

\ m=l 



J_ c| 1 - j f m Z 

N eg 1 - i r m z 



In the following we show that for all z G Co 

c a nsw (z) = v(3t(*)(z) - 1) + > - 



(57) 



which means that (38) holds. Then, according to Corollary 1 the function 
t — > K(x,t) is causal. 

As in the proof of Proposition 1 we prove r]{^.{^){z) — 1) > and 
£Ss{V){z) > 0. 



• Taking into account (24) we define 

„2 / n 



s :- 



TV 5g 



Vm=l 



(58) 



Using this notation, we get 



& m (z) = s 



= s 



(1 + frn T m \z\ 2 ) +1 (T m Z-f m z) 
\l-\T m z\ 2 

7"m {^~m/^~m ~t~ 7" m \z\ -\- ( y 7~ m /f nl -\- T m 
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Because r m jf m > 1, by assumption (56), it follows that for all z e C 

^ m{z))>s Tm . 

Consequently, by using the definition of s, (58), it follows that 



JV 



JV _ 



£ # m (z) > S ^^ = 1. 

\m=l / m=l Tm 



Now, using (54) it follows 



JV 



\ m=l 



> 



SB f £) *„(*)) > 1 



and consequently ni^ii^/ {z)) — 1) > 0. 
We have 



'm \' m 



T m 1 + \z\ 2 + 2r m ?7 

together with the assumption (56), which state that r m /f m > 1, it follows 
that sgn(3(tf r „ l (z))) = sgn(£). According to our assumption, we take that 
complex root, such that the real part of the argument is non-negative 
which together with property (102) in the Appendix implies 



sgn I 9 



JV 



N 



\ m=l 



sgn 3 



Therefore, 



sgn(3(if(z))) =sgn 9 



JV 



\ m=l 



sgn ($t(^2$ m (z)yj =sgn(C) . 



This shows the assertion. 



Our Model 

Theorem 8. For a , r > and 7 e (1,2] let a^ sh be defined as in (28). 
Then K, as defined in (39), is causal. 
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Proof. The function 5^1 + (— ir z) does not vanish in C . Thus the holo- 
morphic extension of w — > a* ksb (uj) is given by 

- ■ r- <n , * <-\ ao(-i^) 



C 0V /l + (-iT z)7-l ■ 

In the following let z e Co- For proving (38) we make a variable transformation 

* /.n _ —ITqZ _ OtQ — 1Z 

akSb[Z ' ~ y/1 + (-iTofy-l ~ T C ^l + (-i^l) ' 

and define 

&(z) = -L and §{z) = 1 + (-i z) 7 " 1 . 

Then, with this notation, in order to prove causality of K, it suffices to prove 
that for all z £ Co 

— a ksb (z) = V $t(*(z)) + &(*(z)) > . (59) 

As in the proof of Proposition 1 we show that both terms rfSt{^{z)) and 
£$*(<?(z)) are non-negative, and then from Corollary 1 the assertion follows. 

In order to prove (59) we note that the function ^ here is the same as in 
(50) in the proof of Proposition 1 when a is set to 1/2. Thus we can already 
rely on the series of manipulations for <P developed in the proof of Proposition 
1. 

• Since r\ > it suffices to show that ^(^(z)) > 0. We note that for a 
complex number a + i b 

K^) = K^) = ^ R(o+i6) - 



(^f<Hz)j > in C . Since 3? (&{z)j > in C for 7 £ (1, 2], it follows 



Taking into account the definition of it therefore suffices to show that 

K V 

that 3? (^f¥(z)j > in C . 
Now, using that 

1 \ /a-i&\ 1 ... 

a + 1 b J \a z + b z J a z + b z 

it suffices to show that -£^{^W{z)) > for proving that £$S(\P(z)) > 0. 
The proof is along the lines as the analogous part in Proposition 1 by 
taking into account that here 7 £ (1,2) (in Proposition 1 7 £ (0, 1)). In 
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this case we have now that sign of S is exactly opposite as in the proof of 
Proposition 1, which in turn gives that has the opposite sign as 

well, and consequently — &{ z )) > 0. Thus the assertion follows from 
Corollary 1. 

In experiments it has been discovered that several biological tissues satisfy 
a frequency power law (16) with exponent 7 <G (1,2) (cf. [WOO, BRBP010]). 
However, as it has been shown in Theorem 3, such models are not causal. Our 
proposed model approximates the frequency power law for small frequencies, 
which is actually the range where it has been experimentally validated. So, 
our proposed model, is valid in the actual range of experimentally measured 
data and extrapolates the measured data in a causal way. Figure 4 shows a 
comparison of a p i and a^sb in an experimental frequency range. 




frequency in MHz frequency in MHz 

Fig. 4. For 7 = 1.5: Comparison w — > aksb(u) (as denned in (28)) where ao := 
2coto/ 1008(1-7) I (dashed line) and the power law a p i(cj) — |tolj| 7 (as defined in 
(16)). For liquids: r = 1CT 6 MHz (left picture) and for gases: r = 1CT 4 MHz (right 
picture) (cf. [KFCSOO]). Experiments for determining the power law coefficient are 
performed in the range — 60 MHz (cf. e.g. [S95]), which is the basis for the range 
of the represented data. 



Model of Greenleaf and Patch 

Proposition 2. For a n > let a* p be defined as in (25) with the specified 
values 7 <G {1,2}. Then K, as defined in (39), is not causal. 

Proof. For the two specified models we have a gpl (u)) = a \uj\ > and 
oi gP 2 (w) = a uj 2 . The respective holomorphic extensions are given by (cf. Proof 
of Theorem 5 and Theorem 3) 

z £ C 1-4- a* 1 (z) := a z + i — zlog [ — ] (w 7^ 0, fixed) 

7T V^o/ 
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and 

z £ C &* gp2 (z) := a (-i zf . 

The assertion for 7 = 1 follows as in the proof of Theorem 5 and the second 
assertion follows from Theorem 3 for 7 = 2. 

Model of Chen and Holm 

Theorem 9. Let < ai < 1, 7 £ (0,2) and Q as in (26). Then there does 
not exist a constant c > such that 

su PP (g(;t))CB ct (0) for t>0, (60) 

i.e. for each c > the function 1 1-> £/(x, i + |x| /c) is not causal. 

Proof. Let t > be fixed. Assume that x i->- C?(x, t) has support in B ct (0) 
for some c > 0. Then according to the Paley- Wiener- Schwartz Theorem 
(Cf. [GW99, H03]) the map k ^ T~ x {£?} (k,i) is infinitely differentiable. 
We show that this is not possible. According to (26) and (27), we have 

with 



A(k) :=-a lCo |kp , B(k) :-coV|k| 2 -"? |k| 27 . 
Since 7 € (0, 2), the function k i-> exp (A(k) t) is not infinitely often differen- 
tiable at k = and since the holomorphic function sm ^^- ) ^ does not vanish 

at k = 0, it follows that k i-> J"^{£}(k, t) is not infinitely often differen- 
tiable at k = 0. Consequently, x Q(x, t) cannot have compact support, 
which concludes the proof. 



5 Integro-Differential Equations Describing Attenuation 

In the following we derive the integro-differential equations for the attenuated 
pressure p a tt for various attenuation laws. Thereby, we first derive equations 
which the according attenuated Green functions Q (cf. (9)) are satisfying, and 
then, by convolution, we derive the equations for p a tt- The integro-differential 
equations are general in the sense, that they apply to arbitrary source terms 
/, and in particular to the source term / (6) of the forward problem of pho- 
toacoustic imaging with attenuated waves. 

For this purpose, we rewrite V 2 J r {^} by using its definition (9), i.e. Q = 
K * t Go, and the product differentiation rule, which gives 



V2tt (61) 
=V 2 T {K} ■ T {Go} + 2V^ {K} • V J" {Go} + ? {K} ■ V 2 T {Go} . 
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To evaluate the expression on the right hand side, we calculate V.F{if} and 
V 2 F{K}. From (10), it follows that 



VT{K} = -&*'■ J" {if} -sgn, 



(62) 



where /?*' denotes the derivative of /3*(r, co) (cf. (10)) with respect to r. This 
together with the formula (103) in the Appendix implies that 

V 2 J-{if} 
= -V - OS*'- J" {if} -sgn) 

= - (V • sgn) • /?*' • J" {if} - (sgn • V/3") • J" {if} - (sgn • V.F{if}) • /?*' 
2 

w 



~ ■ (3*' - p*" + {(3*'f 



J" {if}. 



(63) 

Inserting (62) and (63) into (61) and using again the identity Q = if * t Go, 
shows that 



1 



2 7T 



V 2 T{G} 



/2n 



~ ■ /3*' - /3*" + (/3*'f 



(64) 



- 2/3*' • J" {if } • (sgn • V J" {G a }) + T {if} ■ V 2 F {G a } . 
From this identity, together with the two following properties of Go, 



VT{Go} 



iuj 1 
co |x| 



J 7 {Go} • sgn, 



and 



it follows that 



V 2 T{Go} + — T{Go} 



2tt 



(65) 



(66) 



V 2 T{G} 



^- ■ /3*' - p*" + {p*' f 



- 2 



LO 



iu> 1 
co |x| 



f3*'-T{G} 



(67) 



-T{K}-5 X . 

Inserting the identity T {if} (x, co) ■ S x = J '{if} (0,co) ■ S x in (67) gives the 
Helmholtz equation 
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7 2 



CO 



(3*" ■ T {Q} - T {K} (0,-)-4 
1 



(68) 



exp(-/T(0,u;))-<5 x . 



To reformulate (68) in space-time coordinates, we introduce two convolution 
operators: 

D.f:=K.* t f and LfJ := K' t * t f, (69) 
where the kernels and K' r are given by 

K* :=K*(x,t) :=K,(\x\,t) and K*(r, t) := {/?*'} (r, t) (70) 

V 27T 



and 



^:=^(x,t):=^:(|x|,t) and K' t (r,t) 



T~ x {/?*"} (r , t). (71) 



Using these operators and applying the inverse Fourier transform to (68) gives 



w 2 g 



L21 

co dt 



(72) 



In the case that /3*(|x| ,co) = a*{oj) |x| is of standard form (13), it follows that 
1 



K*{t) = -j=T- 1 {a*}{t) 



and K = 0. (73) 
For a general source term /, we denote the attenuated wave by p a tt- That 



is 



Patt := Patt(x, t)=Q * x>t / =: Af , 

where A is the convolution operator according to the Green function Q. This 
then shows that p att satisfies the integro-differential equation 



2 1 d 2 p &t t . 



(74) 



where A s denotes the space-time convolution operator with kernel 

K s := K,(x, t) := -(Bg)(x, i) + (D^)(x, t) + if (0, i) • <5 x (x) (75) 

and 

B:=Dt + ^D°. (76) 

Co OT 

Equation (74) is called pressure wave equation with attenuation coefficient 
(3* . We emphasize that (3* determines the operators D* and D' t which in turn 
determine the operator A s , which in turn determines p att - this reveals the 
dependence of p att from /3* . 
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Remark 4- Let /?*(r, w) = a*(w)r be the standard attenuation model (cf. (12)). 
Assuming that the associated kernel K (cf. (39)) is causal, it follows that 

\VK\ = I J- 1 {a* • exp (-a* |x|)}| . 



'2tt 

Using some sequence {x„} satisfying x„ ^ and x„ — > shows that 
Jim |VA-| (x n ,i) = -L l^- 1 {a*} (t)| ^1^.(4)1 • 

(73) 



2tt 



Due to the causality of K the left hand side is zero for t < 0, and thus if* is 
also causal. 

Because the convolution of causal distributions is well-defined, the operator 
is well-defined on all causal distributions. Moreover, since = 0, it 
follows that = 0. Using that depends only on t it follows that 

(£>*£) * x ,t / = [K* * t Q\ * Xit / = K* * t [Q * x , t /] = *x,t /)• 

Convolving each term in (72) with a function /, using the previous identity 
and that D' % = 0, it follows that 



V 2 p a tt 



1 <9 n:: 



c at 



Patt - -/ (77) 



where 

I».- = -^=^{a*}(t)* t .. (78) 

In the following we derive the common forms of the wave equation models 
corresponding to the various attenuation models listed in Section 2. 



Power Laws 



• Let < 7 ^ N and < a - We note that the Riemann-Liouville frac- 
tional derivative with respect to time, denote by (see [KST06, P99]), 
is defined in the Fourier domain by 

F{Dlf} = (- iw )^{/} , (79) 

and satisfies 

Dpf = D]D]f and ^D] f = D] ±f = D] +1 /. (80) 



From this together with (15) and (78), we infer 
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and thus wave equation (77) reads as follows 



V 2 p att 



- TVY 19 

c at 



Putt 



(82) 



• Let 7 = 1, then for the frequency power law (17) it follows from the Fourier 
transform table I in [L64] 



D* = — 



^-(logM)Jf 



4ao 

4a H(t) 4a ,. 9 

— * 4 + 7f^ (log|wol) ^- 



Szabo's Attenuation Law: 

Let < a and < 7 £ N. From (19) and (79), we get 



1 a 

Co OT 



1 ^ 2q 

c 2 at 2 c ot * 



and thus wave equation (77) reads as follows 



VVtt 



1 <9 2 p att 2a 9 7 . , 

Cn Ct z c at 



(83) 



Thermo- Viscous Attenuation Law: 

From (20) we get 



Td + T0 | 



1 d 
D* -\ — 

Co ot 



1 2 



liL. 

c 2 at 2 



and thus (77) becomes 



Patt_ ^ 



C5 at 2 



Id + .o-)/. 



This equation is called the thermo-viscous wave equation. 



(84) 



Nachman, Smith and Waag [NSW90]: 

We carry out the details only for one relaxation process. 
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2!) 



N = 1: From (22) we get 



a*(w) + 



(-iw) 



co 



(— i cj) 2 1 — i fi w 
2q 1 — i t\ uj 



which implies 



Id + n 



dt 



c 9i 



Id + n 



d\ 1 d 2 



dt ) c 2 dt 2 ' 



Thus (77) reads as follows 




(85) 



If the term with fi = is dropped and Co is replaced by Co, then we obtain 
the thermo- viscous wave equation (84). 
• N > 1: For the general case we refer to equation (26) in [NSW90]. 



Greenleaf and Patch [PG06]: 

• For 7 = 2 the attenuation coefficient equals to 

a*(u}) — a a u} 2 = 6e (— iuj) 2 
where So is defined as in (81) and thus 

d 2 



D* = -ao 



dt 2 



which gives wave equation (82) with 7 = 2. 
• For 7 = 1 we have 



and thus 



a*(uj) = a \u\ = a (— iuj) isgn(cj) 



d 

D* = a D 1-1 = -ce ttlH , 
at 



where D _1 and H denote the Riesz fractional differentiation operator and 
the Hilbert transform (cf. Appendix), respectively. Therefore the wave 
equation reads as follows 



V 2 p att 



a T>- 1 + — — 

Co at 



Patt 
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Chen and Holm [CH04]: 

Let 7 <G (0,2). The Green function denned by (26) satisfies the Hclmholtz 
equation 



d 2 T{G} 
dt 2 



(k,t)+2a lC0 kf ' V/ ' { - ;| ' > ^ - 2 '^ 



dt 



(k,t) + c 2 |k| 2 T{G}(k,t) 



"(2tt) 3 / : 



■S(t) 



(86) 



for t € M and k € K 3 . Since the fractional Laplacian for a rotational symmetric 
function / and 7 G (0,2) is defined by (cf. Definition (2.10.1) in ([KST06]) 

(-V 2 ) 7/2 /(x) 

:= — ; / cxp(x-k) Ikl 7 — ; / cxp (— x • k) f (x) dx 

V^) 3 J 711 v^) 3 J PV ' ^ ' 

R3 L R3 

=^ 3D {|kr^{/}(k)}(x), 
we obtain the following wave equation for p att := 5 * x ,t / 



dk 



1 <9 2 p att 2ai 9 



(87) 



We note that Chen and Holm used instead of 2ai/co the term 2ai/c 

Our Model [KSB10]: 

From (28) we get 



1-7 



(id+T-r 1 !?- 1 ) 



)_d_ 

Co <9i 



n 2 



1^1 

c 2 di 2 



(aoId + L 1 / 2 )' 



where the time convolution operator L 1 / 2 is the convolution operator with 
kernel 

l(t) := L^iSt) = -jL=r{y/l + (-iro^" 1 } . 

Consequently, L := (L 1 / 2 ) 2 = Id + tq D^ 1 and (77) can be rewritten as 
follows 



2 #2 



Patt 



at 2 

(id + rj- 1 ^- 1 ) /• 



6 Pressure Relation 
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In this section we derive the relation between p att and p when the source term 
is of the form (6). This chapter is a special instance of Section 5. However, 
utilizing the special structure of the source term different formulas can be 
derived. 

Attenuation is defined as a multiplicative law (in the frequency domain) 
relating the amplitudes of an attenuated and an unattenuated wave initial- 
ized by a delta impulse. Here we are concerned in deriving the convolution 
relation between the solution p of (1) (or cquivalently of (4) and (5)) and 
the attenuated wave function p att , which, according to (9) and (29), is given 

by 

Patt = S*x,t/ = (#*t 0O) **,*/, (89) 

with / from (6). Using (7) and the rotational symmetry of K, it follows that 

Patt(x,*) 

(K * t 0o)(x - x', t - t")p(x') d^^{t") dt" 

J J J K(x- x ',t-t' -t")g ( x - x ',t')dt' P ( x ')d x '^(t")dt" 

R R3 R 

J^(x - x', t - t' - t")£o(x - x', t')p{x')6 t (t") dt' dx' dt" 
= /| K (| x - xltH x-x1/c.) i -^L I d x '. 

R3 

Using the representation x' — x = r's with r' > and s e S 2 , it follows 
that 

oo 

p att (x,i) = i- J ^- t K(r',t-r'/c )r' J p{x + r' s) ds dr' . (90) 

o s 2 

Moreover, 



32 Richard Kowar and Otmar Scherzer 



d f 6 t (t- |x-x'|/co) , ,s , , 
«>(*,*) = -+J 47 r|x-x'| 



dt J 47T |x — x' 

oo 

d_ f , 2 fS t (t-r'/c ) 
dt, 



J, 2 jSji^M p{x + r , s)dsdr , 

5 2 

oo 

= di J r 4^ J p(* + rs ) dsdr 

o s 2 

oo 

= ^J ^ 6 t{t-r") J p(x + c r"s)dsdr 



s 2 

dt 4tt 



^ J p(x+ (c Q t)s) ds 
s 2 ) 



This gives 



r' Icq 



S 2 



r' J p(x + r's) ds = J p (x,t')dt' 

Now, denoting 

r' 

F(t,r'):= J ^-K(r",t-r"/c )dr", G(r') := |p (> 



r'/c 

(X. 



and 



F(t,oo) = lim F(t,r') 
it follows from (90) and (91) and the fact that po(-x,0) = that 
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CO 

p at t(x,t) =- / F'(t,r')G(r')dr' 

CO J 



oo 

= - 1 / F(t,/) G'(r') dr' + 1 F(t, r')G(/)|~ =0 
Co J y — v — ' c 



=Po(x,r'/c )/c 

F(i,oo) / p (x, r'/c )dr' - I F(t,r')p (x,r /c ) dr' 



(OO oo > 

F(i,oo) y p (x,t')dt' - J F(t,c t')p {x,t')dt' 


oo 

: y M(t,t')p (x,t')dt' , 



(92) 



where 



t') : = — oo) - F(t, c t')) ■ (93) 
Co 

In the following we derive an equivalent representation of M. in terms of the 
attenuation coefficient, under the assumption that the attenuation coefficient 
a* is such that K is causal. From (10), (12) and Item 4 in the Appendix, it 
follows that 

K(r', t - r'/co) = -j= T~ x jcxp (-a*(w) r' + i ^ r'^j J (t) . 

which implies 

F{t,c t') 

c t' 



J ^-K{r',t-r'/c )dr' 



1 Jzr _ 1 f -i a; [cxp (-a* (a>)c t' + i u;t') - 1] ^ ;i . 



(94) 



^/2tt \ -a*(u>) + i lo/cq 

1 _ i f -i w 

Since X is causal, it satisfies K(cot' ,t — t') = for t < t', and therefore 



•r'{ ^M? lu/a , }w«W'>>'-*-*)-m- 



Fit, oo) = lim Fit, c t') = — L= T~ x \ — —. — 1 it) . 

V ' *'->«> v ' u ; \-a*(cj) +i lo/c J w 

Hence (94) can be written as follows: 
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F(t, c 1') = F(t, oo) * t K(c a t', t-t')- F(t, oo) . 
and therefore (93) simplifies to 

M{t,t') = -—F(t,oo) * t K(c t',t-t') 

Cq 

1 jr _ 1 f iuj C xp(-a*(Lo)cot> + iojt>) \ (f) (95) 



y/2ir \ -a*(oj) c + i uj 

Note that A4(t, 0) = — ^F(t, oo). The following lemma shows that if if is 
causal 

M(M') = if < i < t' (96) 

and therefore the upper limit of integration in the last term (92) can be 
replaced by t. This means that the set of attenuated pressure values 

{p att (x,s):0<s<t} 

depend only on the unattenuated pressure values 

{p (x,s) : < s <t} . 

Lemma 3. Let K from (10) be causal with /3*(r,u>) :— a*{uj)r. Moreover, 
assume that a*(u>) ^ i uj/cq for every and let M. be as defined in (95). 

Then, 

• the function t — > A4(t, 0) is causal. 

• For every t' > 

M(t,t') = for allt < \t'\ . (97) 

Proof. • In order to prove causality of t — > M. (t, 0) we verify the three 
assumptions of Theorem 2 for the tempered distribution 

F{M}(lo,0) = -L= — - s with k(w) :=ia*y + - . 

V2^ k(uj) c c 

Since K is causal, as has been shown in Remark 4, also the function 
t^K*(t) = -L=T- 1 {a*}(t) 



is causal. Now, using Theorem 2, it follows that 

1) a* is holomorphic in Co, 

2) + i rj) a*(£) for n -> in 5' and 

3) for each e > there exists a polynomial P such that |a*(z)| < -P(|z|) 
for z e C e . 
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Since a*(w) ^ i w/c for all w e 1 together with in [BW66, Theorem 
2.7] it follows that z — > a^(z) is unique holomorphic extension to Co and 
therefore z — > a*(z) cannot be identical to z — > i z/cq, the holomorphic 
extension of i w/c . a* (z) ^ i z/c for z e Co implies that k has no 
zeros and hence zjkiz) is holomorphic on Co- This shows that Item 1 in 
Theorem 2 is satisfied for t -> A4(i, 0). 

Since l/(fc(£ + i 77) fc(£)) is bounded and fc(£ + i rf) -> fc(£) for 77 -> in S", 
it follows that for all V e «S 



lim 



/ 



i.e. + i 77) -> l/fc(£) for 77 -> in S". Hence (£ + i 77, 0) -> 

(£, 0) for 77 — > in 5', which shows Item 2 in Theorem 2. 
Since fc(z) does not vanish on Co and |fc(z)| is bounded by a polynomial 
in \z\ for z G Co, it follows that |z/fc(z)| is bounded by a polynomial in 
\z\. Hence Item 3 in Theorem 2 is satisfied and consequently t — > A4(t, 0) 
is causal. 

Property (97) is satisfied if 

M(t+|i'|,i') = for t<0. (98) 



From (95) and 



K(c t',t)=F- 1 \ =— \(t). 



/2ir 



it follows that 



M(t + \t'\,t') =M(t,0) * t K(c t',t). 

Since t i->- A4(t, 0) and i M> K(c n t',t) are causal, their convolution is 
also causal (cf. Item 7 in the Appendix). This proves property (98) and 
concludes the proof. 

Remark 5. Assume that the attenuation coefficient is given by 

a*(u>) = i uj/cq for ui € R . 

Then 

K(r, t) = —= J- 1 {cxp (-i w r/c )} (t) = <J(i + r/c ) 

V 2 7T 

which implies together with (9) and (7) that 
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But this function does not correspond to the intuition of an attenuated wave, 
which is manifested by the convolution equation (9), which should give a 
smooth decay of frequency components over travel distance. With this Green 
function Q the input impulse collapses immediately and consequently, in this 
case, the assumption a*(cj) ^iw/c in Theorem 3 reflects physical reality. 

7 Solution of the Integral Equation 

The inverse problem of photoacoustics with attenuated waves reduces to solv- 
ing the integral equation (92) for p , and to the standard photoacoustical in- 
verse problem, which consists in calculating the initial pressure p in the wave 
equation (4) from measurements of Po( x > t) over time on a manifold surround- 
ing the object of interest. The standard photoacoustical imaging problem is 
not discussed here further, but we focus on the the integral equation (92). 

In the following we investigate the ill-conditionness of the integral equation 
(92), where the kernel M. is given from the attenuation law (28) with 7 <G (1, 2]. 
In this case the model is causal and the parameter range 76 (1,2] is relevant 
for biological imaging. 

In order to estimate the ill-conditionness of the integral equation (92) it 
is rewritten in Fourier domain: 



After discretization the ill-conditionness of this equation is reflected by the 
decay rate of the singular values of the matrix T {M} (w, t') at certain discrete 
frequencies and time instances. 

We consider simple test examples of attenuation coefficients p (as in (3)), 
which are characteristic functions of balls with center at the origin and radiii 
R. For these examples we investigate the dependence of the ill-conditionedncss 
of (99) on the radius R and the location Xo. For applications in photoacoustic 
imaging x would be the location of a detector outside of the object of interest, 
to be imaged. Then, by solving the integral equation (99) po can be calculated, 
and in turn, the absorption energy p can be reconstructed with standard 
backprojection formulas. Since 




(99) 



E 



supp(p (x , •)) = [(l x o| - R)/cq, (|x | + R)/c ] 



the integral equation (99) can be rewritten as 



00 





(100) 



(R -R)/c 



where Rq = |x 
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In the following we analyze the integral equation (100) in terms of the two 
parameters R and R = |x |. This gives a clue on the effect of attenuation in 
terms of the size of the object and the distance of the location x to the simple 
object. In order to show the effect of attenuation on the single frequencies, 
we make a singular value decomposition of the kernel of the integral equation 
(100). 

Example 2. For small frequencies the attenuation law of castor oil, which be- 
haves very similar to biological soft tissue, is approximately a power law with 
exponent 7 = 1.66 and a w 4 • 10~ 2 — ^ IHz) ~, , i.e. 

a p i(w) ss 4 • 10~ 2 • w 1 ' 66 • cm- 1 (w in MHz). 

The sound speed of castor oil is 1490 • ^ at 25 degree Celsius. In units of cm 
and MHz we have 

c w 0.15 • cm • MHz. 
Since (28) approximates the power law (cf. Figure 4) it follows that 

a sin(f (7-I)) 

2 Co To 

and consequently the coefficients of (28) satisfy 

2 c a . 
ro^ 1 ' sin(f (7-I)) 

We note that the relaxation time is r 10~ 4 M 1 Hz for liquids (cf. [KFCS00]). 

For the calculation of the singular value decomposition of the discretized 
kernel of the integral equation (100) we used a frequency range 
lu € [-80,80]M£/z and step size Au = MHz with N = 2 9 . The time 
interval has been set to [0, ^]MHz^ 1 and a step size At — ^ MHz^ 1 was 
used. 

The upper left picture in Fig. 5 visualizes the discretized kernel of the 
integral equation (100) for i?o = R, l - c -, when xo is directly on the surface 
of the object of interest. The upper right picture shows the singular values of 
the discretized kernel in a logarithmic scale. Two properties of the singular 
values become apparent: 

1. For large indices the decay rate is exponential, which can be seen from 
the linear decay in the logarithmic scale. 

2. Secondly, there is a range of indices, where the singular values do not 
decay that rapidly. As a consequence, for solving the integral equation 
this means that the Fourier coefficients of po according to the first block 
of singular values can be determined in a stable manner. 
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For increasing distance L = |x | — R of x to the object the singular values 
of the discretization of the integral equation (100) show a drastically more 
exponentially decay rate for increasing L (see bottom right picture Fig. 5). 
This means that if the object is further away from Xo attenuation is more 
drastically, and solution of the integral equation is more unstable. We analyze 
the dependence of the number of largest singular values from L. For this 
purpose we denote by n cu t the index of the singular value that is about 0.1% 
of the maximal singular value. For the numerical solution of (100) it means 
that if we make a truncated singular valued decomposition with only fi C ut 
singular values, the error amplification can be bounded by a factor 1000. The 
dependence of n cut on L is shown in the lower left picture of Fig. 5. The 
picture reveals that for increasing distance (from about 2cm) only about four 
Fourier modes of po are significant when a maximal error amplification of a 
factor 1000 is required. This reveals that in general the solution of the integral 
equation (100) is significantly ill-posed and worse if the data recording is far 
away from the object. 

Example 3. An analogous numerical example as in Example 2 for the case 
7 = 1.1 is presented in Fig. 6. From the lower left picture of Fig. 6, we see 
that if the distance is about 2 • cm from the boundary of the object, then 17 
singular values are available for the numerical estimation. 

Example 4- An analogous numerical example as in Example 2 for the case 

7 = 2 is presented in Fig. 7. From the lower left picture of Fig. 7, we see that 
if the distance is about 2 • cm from the boundary of the object, then only 4 
singular values are available for the numerical estimation. 

Example 5. An analogous numerical example as in Example 2 for the fre- 
quency power law 

a. (a;) = a p Q l • (-i w) - 66 (a p Q l as in Example 2) 

is presented in Fig. 8. From the lower left picture of this figure, we see that if 
the distance is about 2 • cm from the boundary of the object, then 77 singular 
values are available for the numerical estimation. If the distance is about 4cm, 
then 46 singular values are available. 

Comparing all numerical examples shows that the larger 7 (stronger at- 
tenuation), the more rapidly decrease the singular values. 

8 Appendix: Nomenclature and Elementary Facts 

Sets: Br denotes the open ball with center at and radius R. S n C K™ 

denotes the n-dimensional unit sphere. 
Real and Complex Numbers: C denotes the space of complex numbers, K. 

the space of reals. For a complex number c = a + \ b a — 5i(c), b — 3(c) 

denote the real and imaginary parts, respectively. 
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For a complex number c we denote by |c| the absolute value and by <f> € 
(— 7r,7r] the argument. That is 

c = |c| exp (i</>) . 

As a consequence, when w — r exp (i <j>) then 

w 1 = exp (7 (log(r) + i </>)) . (101) 

Consequently w 1 has absolute value r 1 and the argument is j(f> modulo 
2ir. In this paper all power functions are defined on C\M_. We note that 

«>eC\L and 3£(V^)>0 => %{w) S(y/w) > . (102) 

Differential Operators: V denotes the gradient. V- denotes divergence, and 

V 2 denotes the Laplacian. 
Product: When we write • between two functions, then it means a pointwise 

product, it can be a scaler product or if the functions are vector valued 

an inner product. The product between a function and a number is not 

explicitly stated. 

Composition: The composition of operators A and B is written as AB. 
Special functions: 

• The signum function is defined by 

sgn := sgn(x) := - — - . 



In R 3 it satisfies 

The Heaviside function 
H 

satisfies 



2 

V • sgn = — . (103) 
x 



-*<«>-{! 



for t < 
for t > 



H := 1(1 + sgn). 

The ^-distribution is the derivative of the Heaviside function at and is 
denoted by St := St(t). In our terminology St denotes a one-dimensional 
distribution. Sometimes, if the context is clear, we will omit the sub- 
script at the ^-distributions. 

The three dimensional (5-distribution <5 X is the tensor product of the 
three one-dimensional distributions S Xi , i = 1,2,3. Moreover, 

S x ,t :=<5x,t(x,t) =S x -5 t , (104) 



is a four dimensional distribution in space and time. If we do not add 
a subscript S denotes a one-dimensional ^-distribution. 
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• Xn denotes the characteristic set of Q, i.e., it attains the value 1 in J? 
and is zero else. 

Properties related to functions: supp(g) denote the support of the function 

g, that is the closure of the set of points, where g does not vanish. 
Derivative with respect to radial components: We use the notation 

r := r(x) = |x| , 

and denote the derivative of a function /, which is only dependent on the 
radial component |x|, with respect to r (i.e., with respect to |x|) by •'. 
Let (3 = /3(r), then it is also identified with the function (3 — /?(|x|) and 
therefore 

|x| 

Convolutions: Three different types of convolutions are considered: * t and * w 
denote convolutions with respect to time and frequency, respectively. Let 
/, /, g and g be functions defined on the real line with complex values. 
Then 

/ *t 9 ■= J f(t - t')g{t')dt', f * w g := J f(cu - w')g(w')duj' . 

R R 

* Xi i denotes space-time convolution and is defined as follows: Let /, g be 
functions defined on the Euclidean space R 3 with complex values, then 

/ *x,t 9 ■= I I /(x - x', t - t').9(x', t')d^dt' . 



Fourier transform: For more background we refer to [L64, T48, P62, Y95, 
H03]. All along this paper T {■} denotes the Fourier transformation with 
respect to i, and the inverse Fourier transform J 7-1 {•} is with respect to 
co. In this paper we use the following definitions of the transforms: 

.f{/}m = -^=y cxp (iut)f(t)dt, 

R 

J- 1 {/} (t) = --= J cxp (-iwi) f(oj)doj . 

R 

The Fourier transform and its inverse have the following properties: 
1. 



2. 
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T if ' 9} = -^^{f} * w F{g} and 

V27T 



T{f}-T{ 9 } = ^T{fng}, 



^{f.g} = ^=^{f}n^{g} and 

3. For a e R 

J- {/(t - a)} (u) = cxp (taw) • J 7 {/(t)} (w) 

4. The ^-distribution at a e K satisfies 

S t (t — a) — -^=T~ X {cxp(i au;)} (t) . 
\/2ir 



5. Let / be real and even, odd respectively, then T {/} is real and even, 
imaginary and odd, respectively. 

6. The Fourier transformation of a tempered distribution is a tempered 
distribution. 

7. Let n, T2 € K. If /1 and / 2 are two distributions with support in 
[ti,oo) and [r 2 ,oo), respectively, then f\ * f 2 is well-defined and (cf. 
[H03]) 

supp(/i * / 2 ) C supp(/i) + supp(/ 2 ) C [n + t 2 , 00) . (105) 
The Hilbert transform for L 2 — functions is defined by 

H{f}(t) = U{^ds, 
ttJ t - s 

R 

where £ R f(s)ds denotes the Cauchy principal value of J M f(s)ds. 

A more general definition of the Hilbert transform can be found in [BW66] . 

The Hilbert transform satisfies 

. H{F{f}}(iu) = -iF{ S gnf}(uj), 

. % {% {/}} = -/. 

From the first of these properties the Kramers-Kronig relation can be 
formally derived as follows. Since f(t) is a causal function if and only if / = 
H-f and H= (l+sgn)/2, it follows that T {/} = [T {f}+iU {T {f}}]/2, 
which is equivalent to = iH {J 7 {/}}, i.e. 

K(^{/}) = -9(W{^{/}}) and Z(F{f}) = Vt(H {F {/}}). 
The inverse Laplace transform of / is defined by 
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where 7 is appropriately chosen. 

The inverse Laplace transform satisfies (see e.g. [H91]) 



ZT 1 {h(s - a)} (t) = cxp (at) ZT 1 {h(s)} (t) for all a, t e K (106) 
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Fig. 5. Case: 7 = 1.66 (castor oil). The upper left and right pictures visualize the 
kernel T {Ml} (w, t') and its singular values for L :— Ro — R — Ocm. The lower 
right and left pictures visualize F {M} (u},t') for the detector distances L = • cm, 
L — 2 • cm and L = 4 • cm and the respective indices n cu t for which the singular 
values are about 0.1 per cent of the maximal singular value. 
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Fig. 6. Case: 7 = 1.1. The upper left and right pictures visualize the kernel 
T {M} (uj, t') and its singular values. The lower right and left pictures visualize 
J-{M} (u>, t') for the detector distances L = • cm, L — 2 • cm and L = 4 • cm and 
the respective indices n cu t for which the singular values are about 0.1 per cent of 
the maximal singular value. 
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Fig. 7. Case: 7 = 1.1. The upper left and right pictures visualize the kernel 
T {M} (uj, t') and its singular values. The lower right and left pictures visualize 
J-{M} (u>, t') for the detector distances L = • cm, L — 2 • cm and L = 4 • cm and 
the respective indices n cu t for which the singular values are about 0.1 per cent of 
the maximal singular value. 
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Fig. 8. Case: 7 = 0.66. The upper left and right pictures visualize the kernel 
T {M} (uj, t') and its singular values. The lower right and left pictures visualize 
J-{M} (u>, t') for the detector distances L = • cm, L — 2 • cm and L = 4 • cm and 
the respective indices n cu t for which the singular values are about 0.1 per cent of 
the maximal singular value. 
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£>*, 26 
K„ 26 
Co, 13 
C e , 13 
£>, 10 
M, 33 
5, 10 

Patt, 26 

absorption energy, 2 
attenuation, 4 
attenuation coefficient, 4 
attenuation coefficient, /3*, 4 
attenuation coefficient, standard, 4 
attenuation law, 4 
attenuation law, [KSB10], a£ s() , 7 
attenuation law, Greenleaf & Patch, 

Q*p, 6 

attenuation law, Nachman & Smith & 

Waag, a* naw , 6 
attenuation law, power law, a*;, 5 
attenuation law, Szabo, a* z , 5 
attenuation law, thermo- viscous, a* tv , 6 

convolution kernel, K, 4 



finite speed of propagation, 9 

Function, causal, 7 

Function, rapidly decreasing, 10 

Green function, 8 

Green function, standard wave equation, 
3 

Operator, causal, 8 
Operator, rotation invariant, 8 
Operator, shift invariant, 8 
Operator, strongly causal, 9 

Riemann-Liouville fractional derivative, 
DJ, 27 

source term, photoacoustic model, 3 

tempered distributions, 5', 11 
travel time of a wave front, 8 

wave equation, 2 

wave equation, photoacoustic model, 2 
wave front, 9 
wave front speed, 9 



